#!/opt/local/bin/python
#-*- Encoding: UTF-8 -*-

"""
===============
cross correlate

..codeauthor.. Ralph Kube <ralphkube@googlemail.com>

Compute the cross-correlation C of two timeseries X(t), Y(t), defined as:

C_XY(tau) = E[ (X_t - mu_x)*(Y_t+tau - mu_y) ] / sigma_x sigma_y
C is even: C_XY(tau) = C_XY(-tau)

Where 
mu_t = E[ X_t ]
mu_y = E[ Y_t ]
are the mean of the timeseries X(t), Y(t) and the expectation value is defined as:

E[ ] = Sum_i=1^N t_i p_i

where t_i is the value of the timeseries at t=i.

"""

import numpy as np

def correlate(X, Y, window_length, mode=1):

    """
    Compute the cross-correlation of signal1 and signal2 within window_length time-lags

    rho_{x,y}(tau) = E[((X(t) - E[X]) * (Y(t + tau) - E[Y]))] / (STD[X] STD[Y])

    """
    
    if (X.size != Y.size):
        raise ValueError('signal1 and signal2 have to be the same length')
    if (window_length > X.size):
        raise ValueError('window_length must be smaller than the signal')
    corr = np.zeros(2 * window_length + 1)
   
    mu_X = X.mean()
    std_X = X.std(ddof=1)

    mu_Y = Y.mean()
    std_Y = Y.std(ddof=1)

    # Compute correlation only for positive time lags., tau is the timelag

    # Correlation at tau=0
    corr[window_length] = ((X - mu_X) * (Y - mu_Y)).mean() / (std_X * std_Y)
    for tau in np.arange(1, window_length+1):
        #print '%d/%d' % (tau, window_length)
        #if mode == 1:
            #corr[ window_length + tau ]     = ( (X[:-tau] -X[:-tau].mean())  * (Y[tau:] - Y[tau:].mean())).mean() / ( X[:-tau].std() * Y[tau:].std() )
        corr[window_length + tau] = ((X[:-tau] - mu_X) * (Y[tau:] - mu_Y)).mean() / (std_X * std_Y) 
        #elif mode == 2:
        #    corr[ window_length + tau ]     = ( (X[:-tau] -X[:-tau].mean())  * (Y[tau:] - Y[tau:].mean())).mean() / ( X.std() * Y.std() )
    corr[:window_length] = corr[-1:window_length:-1]
    
    return corr


def fwhm_estimate( corr, dt ):
    n = np.size(corr)

    for n_up in np.arange(n/2+1, n):
        if ( corr[n_up] < 0.5  ):
            break

    for n_down in np.arange( (n-1)/2, 0, -1 ):
        if ( corr[n_down] < 0.5 ):
            break

    return [n_down, n_up], dt*(n_up-n_down)










